home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / lpc / synthesis / inters.c < prev    next >
C/C++ Source or Header  |  1991-12-20  |  4KB  |  133 lines
  1. /* inters.f -- translated by f2c (version of 26 January 1990  18:57:16).
  2.    You must link the resulting object file with the libraries:
  3.         -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #include "../../H/f2c.h"
  7.  
  8. /* Subroutine */ int inters_(bold, bnew, deltb, npols, b)
  9. real *bold, *bnew;
  10. real *deltb;
  11. integer *npols;
  12. real *b;
  13. {
  14.     /* System generated locals */
  15.     integer i_1, i_2;
  16.     doublereal d_1;
  17.  
  18.     /* Local variables */
  19.     static doublereal aold[576] /* was [24][24] */, anew[576]   /* was [24][
  20.             24] */, rold[25], rnew[25];
  21.     static integer ncom1, ncom2;
  22.     static doublereal a[576]    /* was [24][24] */;
  23.     static integer i, j, k;
  24.     static doublereal r[25], delta;
  25.     static integer iminus, npolsm, npolsp;
  26.     static real sum1, sum2;
  27.  
  28.     /* Parameter adjustments */
  29.     --bold;
  30.     --bnew;
  31.     --b;
  32.  
  33.     /* Function Body */
  34.     delta = *deltb;
  35.     ncom1 = *npols + 1;
  36.     i_1 = *npols;
  37.     for (i = 1; i <= i_1; ++i) {
  38.         ncom2 = ncom1 - i;
  39.         aold[*npols + i * 24 - 25] = bold[ncom2];
  40. /* L10: */
  41.         anew[*npols + i * 24 - 25] = bnew[ncom2];
  42.     }
  43.     npolsm = *npols - 1;
  44.     npolsp = *npols + 1;
  45.     i_1 = npolsm;
  46.     for (k = 1; k <= i_1; ++k) {
  47.         i = *npols - k;
  48.         i_2 = i;
  49.         for (j = 1; j <= i_2; ++j) {
  50. /* Computing 2nd power */
  51.             d_1 = aold[i + 1 + (i + 1) * 24 - 25];
  52.             aold[i + j * 24 - 25] = (aold[i + 1 + j * 24 - 25] + aold[i + 1 + 
  53.                     (i + 1) * 24 - 25] * aold[i + 1 + (i + 1 - j) * 24 - 25]) 
  54.                     / ((float)1. - d_1 * d_1);
  55. /* L20: */
  56. /* Computing 2nd power */
  57.             d_1 = anew[i + 1 + (i + 1) * 24 - 25];
  58.             anew[i + j * 24 - 25] = (anew[i + 1 + j * 24 - 25] + anew[i + 1 + 
  59.                     (i + 1) * 24 - 25] * anew[i + 1 + (i + 1 - j) * 24 - 25]) 
  60.                     / ((float)1. - d_1 * d_1);
  61.         }
  62.     }
  63. /* ......COMPUTE AUTOCORRELATION COEFFICIENTS RTEMP(I) WITH R(0) = 1 */
  64.     rold[0] = (float)1.;
  65.     rnew[0] = (float)1.;
  66.     i_2 = *npols;
  67.     for (i = 1; i <= i_2; ++i) {
  68.         sum1 = (float)0.;
  69.         sum2 = (float)0.;
  70.         i_1 = i;
  71.         for (j = 1; j <= i_1; ++j) {
  72.             sum1 += aold[i + j * 24 - 25] * rold[i + 1 - j - 1];
  73. /* L40: */
  74.             sum2 += anew[i + j * 24 - 25] * rnew[i + 1 - j - 1];
  75.         }
  76.         rold[i] = sum1;
  77. /* L30: */
  78.         rnew[i] = sum2;
  79.     }
  80. /* ......COMPUTE R(0) FROM RTEMP(I) */
  81.     i_2 = *npols;
  82.     for (i = 1; i <= i_2; ++i) {
  83.         rold[0] -= aold[*npols + i * 24 - 25] * rold[i];
  84. /* L60: */
  85.         rnew[0] -= anew[*npols + i * 24 - 25] * rnew[i];
  86.     }
  87.     rold[0] = 1. / rold[0];
  88.     rnew[0] = 1. / rnew[0];
  89. /* ......COMPUTE AUTOCORRELATION COEFFICIENTS R(I) WITH CORRECT R(0) */
  90.     i_2 = npolsp;
  91.     for (i = 2; i <= i_2; ++i) {
  92.         rold[i - 1] = rold[0] * rold[i - 1];
  93. /* L50: */
  94.         rnew[i - 1] = rnew[0] * rnew[i - 1];
  95.     }
  96. /* ......INTERPOLATE AUTOCORRELATION COEFFICIENTS */
  97.     i_2 = npolsp;
  98.     for (i = 1; i <= i_2; ++i) {
  99. /* L70: */
  100.         r[i - 1] = delta * (rnew[i - 1] - rold[i - 1]) + rold[i - 1];
  101.     }
  102. /* ......INSERT TEST FOR SINGULARITY OF NEW AUTOCORRELATION MATRIX */
  103. /* ......COMPUTE NEW PREDICTOR COEFFICIENTS */
  104.     a[0] = r[1] / r[0];
  105.     i_2 = *npols;
  106.     for (i = 2; i <= i_2; ++i) {
  107.         sum1 = r[0];
  108.         sum2 = r[i];
  109.         iminus = i - 1;
  110.         i_1 = iminus;
  111.         for (j = 1; j <= i_1; ++j) {
  112.             sum1 -= r[j] * a[i - 1 + j * 24 - 25];
  113. /* L90: */
  114.             sum2 -= r[i - j] * a[i - 1 + j * 24 - 25];
  115.         }
  116.         a[i + i * 24 - 25] = sum2 / sum1;
  117.         i_1 = iminus;
  118.         for (j = 1; j <= i_1; ++j) {
  119. /* L100: */
  120.             a[i + j * 24 - 25] = a[i - 1 + j * 24 - 25] - a[i + i * 24 - 25] *
  121.                      a[i - 1 + (i - j) * 24 - 25];
  122.         }
  123. /* L80: */
  124.     }
  125.     i_2 = *npols;
  126.     for (i = 1; i <= i_2; ++i) {
  127.         ncom2 = ncom1 - i;
  128. /* L110: */
  129.         b[i] = a[*npols + ncom2 * 24 - 25];
  130.     }
  131.     return 0;
  132. } /* inters_ */
  133.